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Abstract 

We present a computation of the next-to-leading order QCD corrections to the production of three Z 
bosons at the LHC. We calculate these corrections using a completely numerical method that combines 
sector decomposition to extract infrared singularities with contour deformation of the Feynman parameter 
integrals to avoid internal loop thresholds. The NLO QCD corrections to pp — ► ZZZ are approximately 
50%, and are badly underestimated by the leading order scale dependence. However, the kinematic depen- 
dence of the corrections is minimal in phase space regions accessible at leading order. 
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I. INTRODUCTION 



The search for and interpretation of new physics at the LHC will require a precise under- 
standing of the Standard Model. Without accurate QCD predictions and reliable error estimates 
for important observables, mistakes in interpreting experimental results may occur. Notable recent 
examples where poor theoretical understanding has hindered the analysis of an experimental result 
are the deviation of the Brookhaven muon g — 2 measurement from the Standard Model prediction, 
the excess of bottom quark production in Run I of the Tevatron, and the discrepancy in the weak 
mixing angle obtained by NuTeV. At the LHC, all analyses require perturbative calculations to at 
least next-to-leading order (NLO) in a s , the QCD coupling constant, in order to make quantita- 
tive predictions that are free from debilitating theoretical uncertainties. This need for higher order 
calculations at the LHC has been summarized in an experimental "NLO wishlist" of processes for 
which QCD corrections are desired |Q]]. 

A cross section at higher orders in perturbation theory consists of two primary components: 
virtual corrections, in which additional loops are added to the Born-level matrix element, and 
real corrections, in which additional partons are radiated. Each contribution is separately infrared 
divergent. They must be combined and summed over degenerate final states to obtain a finite result. 
Initial-state collinear singularities must be absorbed into the definitions of the parton distribution 
functions. 

Well-developed techniques exist for the calculation of real emission corrections at NLO. How- 
ever, the calculation of virtual corrections to processes with many particles in the final state re- 
mains a challenge. For this reason, our primary focus in this paper will be to develop an approach 
to computing virtual corrections to 2 — > 3 scattering processes. We begin by discussing meth- 
ods currently used to perform such calculations. Well-developed techniques exist for calculating 
one-loop virtual corrections to 2 — > 2 and simpler processes. The matrix elements are obtained 
via a standard Feynman diagram calculation. The tensor integrals are reduced to a basis of scalar 
integrals via a reduction algorithm such as Passarino-Veltman [2]. The basis integrals are then 
computed analytically using a Feynman parameter representation, with care being taken to extract 
all infrared singularities that occur in the parametric integration. These integrals are typically per- 
formed with Euclidean kinematics; after the analytic expression is derived, the resulting logarithms 
and polylogarithms are analytically continued to the Minkowski region. 

Several problems arise when this approach is extended to 2 — > 3 and more complicated scat- 
tering processes. The increase in algebraic complexity alone makes the procedure difficult. The 
singularity extraction and analytic continuation to the physical region are typically done on a case- 
by-case basis for each scattering process; the large number of processes for which NLO QCD com- 
putations are desired implies that studying each separately will be an enormously time-consuming 
task. Furthermore, the standard reduction algorithms introduce inverse Gram determinants multi- 
plying the basis integrals; these coefficients vanish when the final-state momenta become linearly 
dependent, and can become arbitrarily small nearby. Although these spurious singularities cancel 
when all basis integrals are combined, it is difficult to establish this analytically, and they usually 
cause serious numerical complications. Careful studies of the boundaries and extrapolations of 
numerical results from safe phase space regions are typically required to obtain stable answers ftj^. 

Because of these complications, the computation of the NLO QCD corrections to 2 — > n, n > 3 
scattering processes is a difficult and intricate task that typically requires one year or more of 
effort for every interaction considered. These difficulties and the importance of these calculations 
to the LHC physics program have stimulated significant effort to develop new approaches for 
perturbative QCD computations. These include new analytic [4] and semi-numerical [5] methods 
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for evaluating loop integrals. Phenomenological results obtained with these methods include the 
NLO QCD corrections to H + 2 jets [@] and ti + jet [7j] at the LHC, and also the matrix elements 
used for W, Z + jets H. 

Ideally, an algorithm for NLO QCD calculations would be highly automated and would handle 
a large number of processes without the need to consider special cases. This would allow a large 
swath of desired corrections to be computed by a single program running in parallel on many ma- 
chines. Any such approach must confront three main issues: spurious phase space singularities that 
appear during the reduction of tensor integrals, the extraction of soft and collinear singularities, and 
the presence of internal thresholds where analytic continuation is required. An approach that ad- 
dresses the first two issues exists, called sector decomposition [Tol Tllll . It permits a completely 



automated, numerical extraction of infrared singularities from loop integrals. The application of 
this approach to an integral results in a Laurent series in e, the parameter of dimensional regu- 
larization, with coefficients that can be numerically integrated over Feynman parameters. Since 
the infrared singularity structure of a loop diagram in parametric space is completely determined 
by its denominator, no Passarino-Veltman reduction of tensor integrals is needed. Consequently, 
inverse Gram determinants never appear. The basic scalar integral which characterizes a diagram 
is identified and sector decomposed. The tensors become polynomials in the Feynman parameters 
after integration over the loop momenta; they can be treated numerically. 

The remaining issue is the internal threshold structure present in loop diagrams. Thresholds 
occur when the internal propagators go on-shell, and a unitarity cut of the diagram leads to a phys- 
ical scattering process. In Feynman parameter space, the denominator vanishes at these points, and 
is regulated only by the — iO prescription for loop integrals. For iV-point functions, this leads to 
denominators with the behavior 1 / (— iO) N ~ 2 at threshold locations. This is completely unsuitable 
for numerical implementation. An approach for handling thresholds in loop diagrams numerically 
was developed in [12l 131. It entails a contour deformation of the Feynman parametric integrals 



off the real axis and into the complex plane to avoid internal thresholds. The integrals are then 
computed numerically. The choice of the deformation for a given diagram is easily automated. 

It appears to us that the combination of sector decomposition and contour deformation provides 
a framework in which numerical calculations of NLO virtual corrections can be fully automated. A 



similar attitude was espoused in [14]. Our goal in this paper is to test this idea on a realistic 2^3 
scattering process at the LHC. We study the NLO QCD corrections to pp — > ZZZ +X, which acts 
as a background to supersymmetric tri-lepton production and appears on the NLO wishlist in [Q]]. 
We find that the combination of these procedures does indeed appear to be a convenient approach 
to NLO QCD computations. 

This paper is organized as follows. In Section|II]we define our notation and present pp — > ZZZ 
at leading order in QCD. In Section|ni]we discuss our computation of the NLO QCD corrections. 
In particular, we present the algorithm we use for the computation of the virtual corrections. In 
Section [IV] we provide numerical results for pp — > ZZZ + X at the LHC. We conclude in Sec- 
tion M 



H. SETUP AND LEADING ORDER PROCESS 

We consider the production of three Z-bosons in proton-proton collisions, 

p(Pi) + p{P 2 ) - Z(p 3 ) + Z(p 4 ) + Z(p 5 ) + X. (1) 
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Within the framework of QCD factorization, the cross section for this process is 

dc=y~W dx 1 dx2ff 1 (x 1 )f'p(x 2 )dai j ^ 3 z+x(xi,X2), (2) 
id Jo 

where the ff 1 are parton distribution functions that describe the probability to find a parton i with 
momentum xPj in the proton pj. The partonic cross sections der^ are computed perturbatively as 
an expansion in the strong coupling constant a s : 

d^ = d^+(^)dag>+0(c#. (3) 



At leading order in this expansion, only the partonic channel q{pi)+q{p2) — > Z (p 3 ) + Z (j» 4 ) + 
Z(p 5 ) contributes. A representative diagram for this process is shown in Fig. (OQ). There are 
six such diagrams, which can be obtained via permutation of the final-state bosons. We neglect 
diagrams containing the exchange of a Higgs boson. For Higgs boson masses below 2Mz, the 
contributions from these diagrams are small. At next-to-leading order, both virtual corrections 
and additional radiative processes occur; we discuss the calculation of these components in later 
sections. After combining the virtual and real corrections, the partonic cross sections in Eq. © 
contain collinear singularities arising from initial-state radiation. These are absorbed into the 
definitions of the parton distribution functions, as discussed in a later section. 

l -Mvm 5 



WW 4 
WW 3 



FIG. 1: Representative Born-level diagram for qq — > ZZZ. 
The expression for the leading order cross section is 

where the factors |, ~, and ~ are from spin- averaging, color-averaging, and identical particles, 
s = X1X2S is the partonic center-of-momentum energy squared, s = 2 Pi ■ P 2 is the total energy 
squared of the proton-proton collision and Vt^ z denotes the final- state phase space. The matrix 
elements |_M^| 2 have the expansion 

|^| 2 = |M 0) i 2 + e|M 0) | 2 , (5) 

where d = 4 — 2e is the space-time dimensionality in dimensional regularization. The matrix ele- 
ments are simple to calculate using standard Feynman diagram techniques. We use a combination 
of the programs QGRAF [15], FORM illol . and MAPLE to obtain them. The required electroweak 
vertex is 



Zqq: HI ~ ~j=^ (</« + 9als) , 

rpq rpq 

9v = ^--Q q s 2 w} g a = --$-■ (6) 
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Here, = 1 — M^/M§ is the sine squared of the electroweak mixing angle, Tf is the weak 
isospin of the quark q, Q q is the electric charge of the quark q in units of the proton charge, and 
Gf is the Fermi constant. Decomposing the matrix elements using the expansion in Eq. ©, the 
leading order cross section takes the form 



da™ = d*® + ed*® . 
The 0(e) term in Eq. © is needed in the computation of the collinear counterterms. 



(7) 



III. NEXT-TO-LEADING ORDER CORRECTIONS 

The 0(a s ) NLO QCD corrections consist of the following components: 

1. the radiative processes qq — > ZZZg, qg — > ZZZq, and qg — > ZZZq; 

2. the collinear counterterms which absorb the initial-state collinear singularities of these ra- 
diative processes into the parton distribution functions; 

3. the one-loop virtual contributions to the leading order partonic process qq — > ZZZ. 

We discuss the calculation of each component in the following sections, and describe in detail the 
method used to compute the virtual corrections. 

A. Real radiation 



We begin with a discussion of the real radiation processes. We present in detail the process 
qq — > ZZZg, and then note the modifications required when qg — > ZZZq and qg — > ZZZq are 
considered. 

Twenty-four Feynman diagrams contribute to q(pi) + q{p2) Z(ps) + Z{p^) + Z(p 5 ) + g(p g ). 
A few representative samples are shown in Fig. ©. Collinear and soft singularities occur when 
the gluon is emitted from an initial fermion line, indicating that the Laurent series for this process 
begins at 1/e 2 . The expression for the cross section is 



da (1) 



1111 



4 9 6 2s 1 3Z+9 



\M 



(1) 



2 dn 



3Z+g- 



The matrix elements are again simple to obtain using standard techniques. 

1 — >^-fU6Ul? g 1 -MWW 5 1 — »— (WW 5 



WW 5 
WW 4 
WW 3 
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WW 3 

mmr- g 
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FIG. 2: Representative diagrams contributing to qq — > ZZZg. 



To discuss the extraction of the infrared singularities, it is convenient to introduce an explicit 
parameterization of the final-state phase space. In terms of the partonic momenta, the phase space 
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takes the form 

x (2 7 r) d d^(p 1 +p 2 -p 3 -p 4 -p 5 -p g ). (9) 

Since there are no singularities associated with the 3Z phase space, we evaluate it directly in four 
dimensions. We now partition the phase space using 

x (2 I r)V< , l(p 1 +p. 2 -p 34! ,-p ! ,) / df! 3z (p 345 ), (10) 

where ^3^(^345) denotes the 3Z phase space with the sum of the three boson momenta giving p 345 
rather than p 1 + p 2 . We evaluate the gluon phase space in the partonic center-of-momentum frame 
by introducing the explicit four-momenta 



p x = ^ (1, 0, 0, 1) , p 2 = ^ (1, 0, 0, -1) , p g = E g (1, 0, s g , c g ) . (11) 

We use the (^-functions to remove as many integrations as possible, and change variables in those 
remaining so that the boundaries are at and 1 . We arrive at the following expression for the phase 
space: 

/ dQsZ+9 = 2(2vr)^r(l - 7) f dX * dX * [Aia-AOPA^ {l - 9z 2 f~ 2e j dn 3Z (p M5 ). 

(12) 

Here, z 2 = Mf / s, and the expressions for the invariant masses in terms of the hypercube variables 
Ai and A 5 are 

s 3 45 = (P3 + P4+p 5 ) 2 = (l-9z 2 )(l-A 5 ) + 9z 2 , 
sig = (p 1 - Pg ) 2 = -\ 5 (l-\ 1 )(l-9z 2 ), 

S2g = ( P2 - P g) 2 = -\ 5 \ 1 (1-9Z 2 ). (13) 

In writing these expressions we have set the overall energy scale s = 1; it can be restored at the 
end using dimensional analysis. 

The singular terms in the matrix elements come from the following three sources. 

1 . Interferences between diagrams where the gluon is emitted from the quark line with momen- 
tum P i. When the denominator of the off-shell quark propagator s lg in Eq. (fT3l) is combined 
with the phase space in Eq. (fT2l) . the singular structure A^ 1-2 ^! — Ai) _1_e is obtained. 

2. Interferences between diagrams where the gluon is emitted from the anti-quark line with 
momentum p 2 : these lead to the singular structure A5 Ajf 1_e . 



3. Interferences containing the denominator 1 / s\ g / s 2g : if care is taken to sum over only phys- 



ical gluon polarizations in the final state, these contain only the soft singularity A 5 1 ' 
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To extract the singularities as a Laurent series in e, we use the following standard plus distribution 
expansion: 

"ln n A" 



A 



-l+e 



n=0 



where the plus distributions are defined via 

7(A) 



e 
n! 



A 



(14) 



dX 



X 



<?(A) 



1 dX ^-(g(X)-g(0)). 



(15) 



After using these expansions, the cross section for qq — > ZZZg takes the form 



da 



(i) 

qq^ZZ+g 



(16) 



The Ai are integrable, e-independent quantities that contain the complete kinematic information 
of the final state. The 1/e 2 singularities, where A5 — > and Ai — > or 1, cancel against the virtual 
contributions to qq — > ZZZ. The 1/e terms where Ai — > or 1 are removed by the collinear 
counterterms discussed in the next section. The 1/e singularities where A 5 — > cancel against a 
combination of the virtual corrections and collinear counterterms. 

The matrix elements for the remaining real radiation processes qg — ► ZZZq and gg — ► ZZZq 
are identical. They each consist of twenty-four diagrams. The cross section is 



a<T qg^3Z+q 



111 



4(1 -e) 24 6 2s'- M3Z+9 



(i) 



2 dQ 



3Z+q, 



(17) 



where we have used the fact that the gluon has 2(1 — e) physical polarizations in d = 4 — 2e 
dimensions. We extract singularities using the same phase space parameterization and expansion 
in plus distributions discussed above. For this process, only collinear singularities where Ai — ► 
or 1 occur. These are removed by the collinear counterterms described in the next section. 



B. Collinear counterterms 

The radiative processes discussed in the previous section contain collinear singularities that 
must be absorbed into the parton distribution functions. To do so, we begin by expressing the bare 
distribution functions and cross sections in Eq. © in terms of the renormalized ones, 



do- = ^ / dx 1 dx 2 ff 1 {x 1 )ff 2 {x 2 )da ij (x 1 , 



x 2 



(18) 



The renormalized parton distribution functions are related to the bare ones fa via 

fi = ® fj- 

We have introduced the convolution integral 



{f®g)(x)= I dydzf(y)g(z)S(x-yz), 
'0 



(19) 



(20) 
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and we have implicitly summed over repeated parton indices. The functions r^- are given by 



, , pW( x ) 

r i3 (x) = 6^6(1 -x)- ^^L±l. (21) 

7T 



The DGLAP kernels P^ in the MS scheme can be found in |170 ; those required here are 

^W = ?{^(l-x) + lr ^ ] --(l + ,) 

= j{* 2 + (1-*) 2 }. (22) 

To proceed, we substitute Eq. ((21]) into Eq. (TT8l) . equate this to Eq. (O, and solve for the renor- 
malized cross sections d&ij. We expand the renormalized cross sections in the strong coupling 
constant 

d*« = d*g> + (^) d*« (23) 

and employ Eq. © to obtain the following relations between the bare and renormalized cross 
sections at each order in a,: 



da^( Xl ,x 2 ) = da^(x u x 2 )- 

1 f 1 r i 

da^(x u x 2 ) = dcrJJ- } + - / dyP$(y) da$ (x 1 ,x 2 y) + da$ (x 1 y,x 2 ) 

e Jo 1 J 

dag\x u x 2 ) = d<7« + - [ dyP^(y)da^{x 1 ,x 2 y). (24) 

e Jo 

The collinear counterterms that must be added to the perturbatively computed cross sections are 
the integrals in Eq. (l24l) . These are straightforward to compute as Laurent series in e, as we do for 
the real radiation cross sections in Eq. (fT6l) . We note that the convolution variable y maps onto the 
invariant mass S345 in Eq. (fT3l . This makes it simple to check analytically that the singularities in 
the real radiation cross section that occur as Ai — > 0, 1 cancel. We also note that the 0(e) term 
in the leading order cross section in Eq. © contributes to the collinear counterterms at the finite 
level. 



C. Virtual corrections 

Finally, the virtual corrections to the partonic process qq — > VVV must be computed. Forty- 
eight one-loop diagrams contribute to this process; these must be interfered with the six tree-level 
diagrams. A representative sample of virtual diagrams is given in Fig. ©. The expressions for 
the interferences are simple to obtain with standard techniques. We regulate all singularities using 
dimensional regularization, in which scale-less integrals are set to zero. Consequently, there are no 
contributions from self-energy insertions on the external legs. All singularities that remain when 
the one-loop diagrams are combined are infrared in origin. 

The treatment of the virtual corrections to 2 — > 3 processes is the main point of this paper, so we 
discuss this here in detail. We will use the pentagon appearing as the rightmost diagram in Fig. © 
to demonstrate our technique. We begin by neglecting all numerator algebra that appears from the 
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FIG. 3: Representative one-loop diagrams contributing to qq — > ZZZ. 



Feynman rules and Dirac traces when this diagram is interfered with the tree-level diagrams, and 
focus on the scalar topology. The basic scalar integral for this topology is 



1 



d d k 1 



1 



1 



1 



(2ir) d k 2 (k + pi) 2 (k + p 3 + p 4 - p 2 ) 2 (k + p 3 - p 2 ) 2 (k - p 2 ) 



(25) 



The +z0 prescription associated with each denominator has been suppressed for notational ease. 
After introducing a standard Feynman parameter representation, this integral becomes 



1 



?r(3 + e) 

(47^/2 



1 5 

Y[d Xi 5 [ l- 



X j 



[A - zO] 



-3-e 



(26) 



n=l 



where 
A = 



2x 2 (x 3 + Xi) pi-p 3 + 2x 2 x 3 pi • p 4 + 2(xi + x 2 ){x 3 + x 4 ) p 2 -p 3 + 2x 3 (x 1 + x 2 ) p 2 ■ p 4 
-2x 2 (x 3 + x 4 + x 5 ) pi-p 2 - 2x 3 (xi +x 2 + x 5 )p 3 ■ p 4 - x 3 (xi + x 2 + x 4 + x 5 ) M| 



-(x 3 + x±)(xi + x 2 + x 5 ) M%. 



(27) 



We must discuss two features of this integral: the extraction of infrared singular terms and the 
treatment of internal thresholds. 

We begin by considering the singular structure. This integral exhibits infrared singularities 
as various combinations of Xi approach or 1 . A convenient, easily automated prescription for 
extracting such singularities from loop integrals was presented in 111 ill . We summarize here the 
salient features of this technique. 

1. To remove singularities that occur as x\ — > 1, we first split the integral into primary sectors. 
There is a primary sector for each Feynman parameter; for example, the primary sector 
associated with x\ is obtained by making the following variable changes in Eq. (|26l) : 



x j 



2,3,4,5. 



(28) 



The 5-function is then used to remove the integration over x\. All singularities are mapped 
to Xi = by this split. 

2. After forming the primary sectors, all singular terms arise when one or multiple Xi vanish. 
We use sector decomposition to handle the cases where several xi go to zero. We illustrate 
this technique on the following simple example: 



dxdy 



1 

(x + y) 



2+e' 



(29) 
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We split this integral into two regions, X = X 1 + X 2 . The first region has x > y, while the 
second has y > x. In the first region we make the variable change y = y'x, while in the 
second we use x = x y. The integrals become 

Xi = f dxdy X i 2 = / dxdy V ■ (30) 

all singularities now occur only when a single x- t vanishes. 

3. The singularities arising from x^ 1 ^ 6 can be extracted using the plus distribution expansion 
in Eq. (fl~4l) . This yields a Laurent series in e whose coefficients can be integrated either 
analytically or numerically. In some cases it is convenient to modify the Feynman parame- 
terization in Eq. ([27]) to reduce the number of sector decompositions required. 



The Feynman denominator A in Eq. (1271) can vanish in the interior of the Xi integration region, 
as it is clear from the presence of terms with both plus and minus signs. This occurs when the 
internal loop particles go on-shell, and signals the onset of an imaginary part in the integral. The 
— iO prescription regulates these internal thresholds. However, this prescription is not suitable 
for a numerical treatment of the integral. A method that allows internal thresholds to be handled 



completely numerically was developed in II 121 . 11311 . The idea is to deform the contour for the 



Feynman parameter integrations away from the real axis in the direction indicated by the — zO term. 
If the contour is sufficiently far from where A vanishes, then the integration can be performed 
numerically in the complex plane. 

We discuss the application of this technique to the pentagon integral in Eq. (|26t . After splitting 
the integral into primary sectors and sector decomposing the integrand, the integrand denominator 
in each sector is given by a product of Feynman parameters factored out during the sector decom- 
position procedure, and a function A that depends upon kinematic variables. This function may 
vanish in the interior of the integration region and, for this reason, its properties determine the 
desired contour deformation. The function A takes on the generic form 

A = Z + ^ YiXi + \ Xi 3 XiX 3 + \ W ijkXiXjX k + ... . (31) 
i i,j 

The tensors X, Y, and W consist of kinematic invariants and are independent of the Feynman 
parameters. The ellipsis denotes terms of quartic and higher order in the Feynman parameters x«. 
Terms beyond quadratic order appear only for 5- or higher-point integrals, and only after the sector 
decomposition is performed. We find that it is necessary to perform the sector decomposition 
before deforming the integration contour; reversing the order can lead to thresholds regulated only 
by — iO. The idea is to now set Xj = yi — iTi, and choose such that internal thresholds are avoided 
and the integration over y { can be done numerically. A convenient choice for Tj, similar to that 
presented in [13], is 



r, = Ajfc(l -jji) 



j j,k 



(32) 



The endpoints of the contour remain fixed with this choice. The parameter A controls the overall 
size of the deformation. Because A contains cubic and higher polynomials in x^, the deformation 
in Eq. (|32l ) does not ensure a sign-definite imaginary part. However, as A — » the deformation is 
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in the direction required by the — zO prescription; this allows us to begin with a small choice of A, 
and check that no pole in the complex plane is crossed as we increase A. For numerical purposes, 
it is convenient to choose A large as compared to the kinematic invariants in X, Y, and W. It 
is simple to obtain the kinematic matrices using computer algebra techniques, indicating that the 
contour deformation procedure can be completely automated. 

Computing the virtual corrections is simple once the singularity and threshold structures of the 
base scalar integrals are regulated. The numerator algebra which arises from computing a com- 
plete diagram has the generic form J\f(k 2 , k- pi) in terms of the loop momentum. After analytically 
integrating over k using standard techniques, the numerator becomes a polynomial in the Feynman 
parameters; it can be treated numerically. No reduction of tensor structures is needed. We found 
it convenient to put each diagram over a common Feynman denominator; large cancellations be- 
tween tensor and scalar integrals occur if the tensor integrals are sector decomposed separately. 
Each interference between one-loop and tree diagrams becomes a Laurent series in e, with coeffi- 
cients that can be integrated numerically. Judicious grouping of terms allows the expression size 
for each diagram to be kept relatively compact. 

In summary, our procedure for computing the virtual corrections is as follows. 

1 . Compute the interference between a one-loop diagram and the tree diagrams using standard 
techniques. 

2. Identify the base scalar integral for each diagram. Introduce a Feynman parameterization 
for this integral, and combine all tensor structures into a numerator over the denominator 
of the scalar integral. Perform the integration over loop momentum analytically; the tensor 
terms become polynomials in Feynman parameters. 

3. Split the base integral into primary sectors, and sector decompose each primary sector until 
all singularities are extracted. The resulting expression will be a Laurent series in e with 
integrable coefficients. 

4. Deform the contour in each sector by making the variable change in Eq. (l32l) . This regulates 
all internal thresholds and allows for a numerical evaluation of the integral. 

IV. RESULTS 

In this Section we describe the results of our computation of the NLO QCD corrections to the 
tri-boson production process pp — > ZZZ at the LHC. We employ the following numerical values 
for the Fermi constant and the weak boson masses: 

G F = 1.166 x 1(T 5 GeV" 2 , M w = 80.451 GeV, M z = 91.1875 GeV. (33) 

We use the MRST lfl8ll parton distribution functions at either LO or NLO, as appropriate. The val- 
ues of the strong coupling constant a.JMz) appropriate to use with the MRST parton distribution 
functions are also obtained from Ref. ft 1 311 . 

We compute the real emission corrections using the procedure described in SectionUm To com- 
pute the virtual corrections, we generate a set of ten thousand random kinematic events distributed 
according to the Born-level matrix element and the leading order parton distribution functions. 
Each event is described by seven variables that provide a complete description of the qq — > ZZZ 
kinematics. We compute the NLO virtual correction for each event and then reweight the event 
using the computed correction and the ratio of NLO to LO parton distribution functions. 
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For the numerical computation of the NLO corrections, we employ the adaptive Monte Carlo 
integration program VEGAS as implemented in the CUBA library [ 1911 . The numerical stability 
of the computation is exceptional; all diagrams including the pentagons exhibit a very fast rate of 
convergence. The computation of the NLO QCD corrections for ten thousand kinematic points 
required a few days of running on a cluster of several dozen processors. 

As an example of our code output we present below in TableUa listing of the finite contributions 
coming from the one-loop corrections for several sample events. The kinematics are defined by 
xi,X2, the values of Bjorken-x for each proton, and the kinematic invariants = (pi — Pj) 2 - The 
initial partonic momenta are denoted by pi,p 2 , while P3,Pi indicate the final-state Z-momenta. 
The invariant masses have been scaled by l/{x\X 2 s) so that their magnitude is between zero and 
one. All other invariant masses can be obtained via momentum conservation. The quantity z = 
Mz/ \jx1X2s is included for completeness. Each event has unit weight at leading order. The shift 
in the weight coming from the NLO virtual corrections is obtained via 

PDF(NLO) / a ST \ 
w NLO = „^, T ' ( 1 + —V) . (34) 



PDF(LO) V vr 

V is the result from numerically evaluating the loop corrections, and is given in the table for each 
event. These must be combined with the real corrections and collinear counterterms to produce the 
final result. We generate these remaining contributions in our program as separate events. Both 
the parton distribution functions and a s can be evaluated at the desired scales. We note that the 
events have been generated using the factorization scale choice jjt,p — 3Mz at leading order. 



{xi,X 2 , Z, S13, S23, S14, S24} 


V 


{0.182, 0.030, 0.088, -0.640, -0.143, -0.132, -0.802} 


7.59(2) 


{0.138, 0.006, 0.226, -0.293, -0.240, -0.487, -0.095} 


9.76(2) 


{0.024, 0.193, 0.096, -0.035, -0.569, -0.692, -0.065} 


9.57(7) 


{0.032, 0.052, 0.160, -0.423, -0.281, -0.226, -0.152} 


7.93(2) 


{0.014, 0.074, 0.202, -0.387, -0.060, -0.040, -0.750} 


11.77(4) 



TABLE I: Finite contributions coming from the one-loop corrections for several sample events. The first 
column indicates the kinematics of the events, while the second column gives the NLO virtual correction 
with its associated Monte Carlo integration error. The notation is as defined in the text. 



We have applied a number of checks to our calculation. 

1 . We have compared the leading order cross-section obtained with our code with the result of a 
similar computation using the program MadEvent [20] and have found complete agreement. 

2. As we mentioned earlier, the NLO virtual corrections are divergent, and physical results are 
only obtained once real emission contributions are added. The divergent part of the NLO 
virtual correction to q(p±} + q(p 2 ) —> ZZZ is related to the leading order cross-section by 
the following equation Bill : 

NLO.virtl n a s r(l + e) [1 3 \ LO 



"""""I- = - C ^1^ (S12) u + fcj ° <35) 

where s 12 = 2p 1 ■ p 2 and Cp = 4/3 is the QCD color factor. We have checked that our 
numerical computation of cr NLO > virt gives the divergent part in full agreement with Eq. (l35l) . 
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3. We have checked that all divergences cancel at the differential level once the real emission 
processes, the collinear counterterms and the virtual corrections are combined. 



4. An important check of the result is provided by its independence of A, the size of the contour 
deformation. However, we stress that the efficiency of the numerical integration depends 
strongly on A. For small values of A one does not move sufficiently far from the pole on the 
real axis, while for large values of A one deforms too much and there are large cancellations 
between different segments of the integration path in the complex plane. 

5. Finally, we have implemented all parts of the computation in at least two independent com- 
puter codes that agree for all observables studied. We have computed the real emission 
processes using both the approach described in the text and the traditional phase-space slic- 
ing method, and have found complete agreement. 



20.0 
17.5 
15.0 
£j 12.5 

b 

10.0 
7.5 
5.0 

2 4 6 8 10 

/x/M z 

FIG. 4: The scale-dependence of the leading order and next-to-leading order cross-sections a(pp — > ZZZ). 
We have set the factorization and the renormalization scales equal to a common value /x. 




We now describe the results of our computation of the NLO QCD corrections to pp — > ZZZ. 
In Fig. |4] we show the dependence of the total cross-section computed through leading and next-to- 
leading order on the renormalization and factorization scales. We have equated these to a common 
scale /x — iir — [Lf- There are two important features of this result to note. The first is that 
the corrections are large, approximately 50% over a wide range of /x. This results from a large 
increase in the qq luminosity function when going from LO to NLO, and large virtual corrections. 
For example, for /x = 3M Z , the LO cross section evaluated with LO parton distribution functions 
is 10.3 fb, while the LO cross section evaluated with NLO distribution functions is 11.4 fb. The 
full NLO result is 15.2 fb, with the additional increase coming entirely from the virtual corrections. 
The effect of real parton emission in the qg and qg channels is 1% or less for all /x considered. We 
note that similarly large corrections for the process pp — > ZZ at the LHC were observed in Il22ll . 

The second important feature is the tiny scale dependence of the LO result, which drastically 
underestimates the NLO correction. The LO result varies by only a few percent over the entire 
range of /x considered. While such behavior is uncommon, it is by no means unique to this process; 
a very similar situation occurs for Z production at the Tevatron Q23D . 



13 




FIG. 5: The transverse momentum and rapidity distributions of the Z bosons at LO and NLO in a s , nor- 
malized by a factor 1/3. The results obtained by re-scaling the LO distribution by a constant A'-factor are 
also shown. The value of the factorization and the renormalization scales are set equal to 3Mz- 



In Fig. [5] we present the transverse momentum and rapidity distributions of the Z bosons. We 
include all three bosons and divide by a factor of 3 to normalize the result. We compare these 
distributions to the approximation of reweighting the LO results by a constant i\T-factor, where K 
is the ratio of NLO to LO inclusive cross sections. For the distributions studied, the NLO QCD 
corrections do not depend significantly on the kinematics of the produced particles. Rescaling the 
leading order kinematic distributions by a constant /^-factor gives a description of the NLO result 
accurate to a few percent. We expect that this is true in all kinematic regions for which phase space 
is available at leading order. 



V. CONCLUSIONS 



In this paper we present a novel numerical method for perturbative computations in QCD to 
next-to-leading order accuracy. We combine sector decomposition lllin . which allows an automatic 
extraction of soft and collinear singularities from virtual diagrams, with contour deformation of 



the Feynman parametric integrals 11121 11311. which permits numerical evaluation of loop corrections 
with internal thresholds. Doing so, we obtain a tool that enables an efficient and flexible numerical 
evaluation of Feynman diagrams with an arbitrary singularity structure. It appears to us that this 
combination of sector decomposition and contour deformation provides a framework in which 
numerical calculations of NLO virtual corrections can be fully automated. 

To test this idea, we compute the next-to-leading order QCD corrections to the production 
of three Z bosons in proton-proton collisions; this is one of the processes that appears on the 
so-called "NLO wishlist" yj]. We observe that the method possesses excellent efficiency and 
numerical stability; for all phase- space points considered, we were able to compute the NLO 
virtual corrections with sub-percent precision. 

The NLO QCD corrections to pp — > ZZZ are large, approximately 50% for all scale choices 
considered. The leading order scale dependence drastically underestimates the size of these cor- 
rections. For phase space regions accessible at leading order, the NLO corrections are independent 
of the kinematics of the final state particles. 
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